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Liquid-liquid-solid transitions (LLST) are known to occur in confined liquids, exist in supercooled 
liquids and emerge in liquids driven from equilibrium. Molecular dynamics (MD) simulations claim 
many successes in forecasting the phenomena. The transitions are also studied in the framework of 
thermodynamics based methods and minimalistic models. In here, the proposed approach is derived in 
the framework of continuum and includes spatial and temporal dynamic heterogeneities; the approach 
is meant to capture the material behavior at small scales. We conjecture that the liquid-like and 
solid-like behaviors are dissimilar enough for the two to be governed by different constitutive relations. 
In this way, we gain additional degree of freedom, which is found essential when predicting the 
transitional phenomena. As a result, we derive the LLST criteria for liquids in equilibrium, during 
steady flow and at transient conditions. Lastly, we forecast short-lived LLSTs in human blood during 
cardiac cycle. 
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Envision a thin film of a molecular liquid exhibiting one or a combination of transitions: the liquid-liquid 
transition (LLT), where a semi-organized molecular structure is formed (liquid II); the liquid-solid trans- 
ition (LST), in this case the liquid is transformed into a solid-like material and, lastly, the liquid-liquid-solid 
transition (LLST) defined as a combination of the two. When controlled, these transitions can be utilized in 
micro-electro-mechanical systems (MEMS), joint lubrication in biology and in micro-fluidics. Suspensions 
exhibit similar behaviors and the transitions play a role in paints, inks, cosmetics, pharmaceuticals and food 1 . 
There are also dense suspensions such as corn starch which exhibit jamming transitions and, in some cases, could 
allow you to run on their surface without sinking. 

The question is: What happens during these transitions? Liquids brought to the transitional regime display 
collective (temporarily and spatially synchronized) motion of molecules and particles. Consequently, relaxation 
times are much longer from those in bulk 24 , viscosity increases by many orders of magnitude and, when sheared, 
confined liquids may experience smooth, stick-slip or chaotic responses 5 9 . The collective motion triggers 
dynamic heterogeneities with a supermolecular length ranging from a few up to ten molecular distances 10-13 . 
Heterogeneities of a similar kind are observed in various suspensions 14 and, among them, in blood subjected to a 
transient flow 15 . We also know that strong stimuli drive viscoelastic liquids to the transitional state. An impact- 
activated solidification has been observed in dense suspensions 16 and an aligned motion of molecules is detected 
in liquids subjected to shock 17 . 

Molecular dynamics simulations provide excellent insight into the transitional phenomena 18-21 . The simula- 
tions are followed by thermodynamics-based models 22,23 and minimalistic models 24 26 ; the latter two are suitable 
for the replication of the observed behaviors. Also, the transitions are studied in the framework of non-Newtonian 
fluid dynamics 27,28 . As we noted earlier, the transitions emerge in over-constrained liquids where molecules lose 
their ability to move freely and, consequently, are forced to act in a collective manner. 

In here, the transitional regime is the regime of our interest. 
Our assertion is that a continuum-level approach can be useful in forecasting the transitions so long as the 
approach is brought close enough to the molecular scale. In our case, this is done by accounting for the relevant 
spatial and temporal fluctuations known as the dynamic heterogeneities 29 . Also, we assume that the substance in 
the liquid-like and solid-like states is different enough for the behaviors to be governed by independent con- 
stitutive equations. By comparing the two we derive criteria for the liquid-solid and liquid-liquid transitions in 
equilibrium and at steady state. The analysis prepares us for a more challenging task, namely the prediction of the 
LLSTs during transient flow processes. We illustrate the later in the example of blood subjected to an idealized 
cardiac cycle. 
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Path Imperfections 




Figure 1 | Trajectory of a material point: starting from the position {X k } 
the new optimal position would be in the point {zjj, but instead the 
particle travels to the actual position {x k }. 

Results 

Dynamic heterogeneity. Let's consider a viscoelastic liquid, where 
the monitored particle moves from its initial position {X k } to a less 
than optimal position {x k }, Fig. 1. The "optimal" position refers to the 
thermodynamically most favorable (mean) positin {z k }. From this 
point of view, a collective motion of a few atoms in mona- 
tomic liquids (transit) can be considered a deviation from the 
expected trajectory 30 . In liquids, the current position of a molecule 
often diverges from the optimal position, but the trajectory must be 
physically admissible, i.e. it must be consistent with the constraints 
imposed by the conservation laws. In the framework of continuum, 
we expect to find the material particle (particle and its surroundings) 
in an acceptable position determined by the equations of motion 
day/ dxj = pbi, where the components of the current stress are Cy, 
the particle acceleration is i), = dVi/ dt and, as usual, mass density is p. 
Any deviation from the thermodynamically optimal trajectory 
[{X k ^>x k } versus {Xk — >z k }] triggers perturbations in stress. 
Consider stress tractions plotted on a surface normal to the 
direction of flow {n k }. The tractions in the optimal and actual 
positions may not be the same a z ik n k i^ a i k n k and the difference is 
responsible for stress fluctuations 



(5V„ 



k )n k njdS. 



(1) 



Stress in the reference (optimal) position is a z ik and in the actual 
position {x k } is a ik . Using divergence theorem, the stress 

perturbations centered about a z k become 5c/-— — UjkknjdV. 

V 0 ) Vo 

The material length l c is understood as the dominant length and it 
captures the relevant spatial stretch of the stress gradient. Next, the 
stress gradient is replaced by the inertia term taken from the 
equations of motion, while volume V 0 is reduced to a material 
point (Vo— >0). As a result, the perturbations are simplified to 
dc/y = l c Oik,k n j — plcYij- We assume that the tensor Sc/y is symmetric 
and, then, we have y« = (o;rtj + ty-M,) /2. However, the symmetry 
restriction does not need to be enforced. The stress fluctuations are 
incorporated into the constitutive description of the liquid and the 
solid. In the transitional regime, the liquid is prone to shear and may 
experience changes in mass density. This behavior is described by the 
Maxwell-like viscoelastic model 



where strain rate is e,j = (pVi/dxj + dvj/dxj) /2. In this relation, elastic 
matrix is Q»j and 1] is viscosity. Viscosity is determined by averages 
over a spectrum of relaxation times 31 . Usually, the viscous term in (2) 
is based on stress deviator alone. In this case, shear stress and 
pressure are viscous quantities 32,33 . Also, the elastic matrix includes 
contributions of bulk and shear moduli. The last term in (2) captures 
the contribution of the stress perturbations. The perturbations <5er^ 
are added to the mean stress denoted as Oy, In the constitutive 
relation, the parameter Rf = V s l c /v k resembles the Reynolds 
number and, for this reason, we call it the material Reynolds 
number. We introduce the number for reasons discussed later. In 
here, kinematic viscosity is Uk = r\j p and u s is sound velocity. When 
the tensors <5o^ and are not symmetric, we expect the non- 
symmetry would trigger perturbations in flow. The rate of 
mechanical work performed by the material is ffi/U;,; = Gj, + 
2\l/ L + R^ctijnjVi/v s . The state function Gi and the dissipation 
potential \j/ L are G L = <T j jC jj ~ k j(T k i/2 and 1^ = 0^^/(2^), 
respectively. We omit the contribution of heat flux. The flux 
R^fy/vs in (2) may become a powerless quantity when the term 
OijYij is equal to zero. As suggested 34 , the powerless flux captures 
the contribution of hidden micro-scale dynamic events. 

The liquid is said to be converted to a viscoelastic solid-like mater- 
ial. In the simplest circumstance, the solid follows the Kelvin-Voigt 
behavior and includes the contribution of the micro-inertia 
described in (1), thus 



Vi 



(s)ds. 



(3) 



The time span (t — t 0 ) is taken to include the relevant history of the 
perturbations. This means that the material retains a short memory 
of the past history but this memory fades away beyond (t — t 0 ) 31,35 . In 



here, mechanical work is Oipij = Fs + 2ij/ s + pv s Vi^rij 



s)ds/R* 



The state function is Fs = Sijdj k i£ k i/2 and the dissipation potential 
becomes \j/ s = ife,- / £^/2. With the use of normality rules 36 the dissipa- 
tion potential captures viscous stress r\&y in (3). The contribution of 
the dynamic heterogeneity (oijtij) p[u;(f) — Vi(t 0 )]v s /R^ is linked to 
the relevant change in particle momentum. Under certain conditions 
the responses produced by the liquid (Eqn. 2) and the solid (Eqn. 3) 
become indistinguishable. This is what we call the liquid-liquid 
transitional state (liquid II). 

LST- near-equilibrium scenario. We place the liquid (2) into a small 
container. Walls of the container restrict the motion of molecules 
and, in this manner, contribute to the increase of viscosity. In 
equilibrium, the fluctuations represent the primary response of the 
substance and, therefore, we omit the inertia terms in (2) and (3). In a 
one-dimensional setting, the nano-scale stress in the position {x} is 
calculated from the equation of motion <t = <j z + p(x — z)b, where a 
and a z are stresses in the actual and optimal positions. At a larger 



(meso) scale, the stress becomes o = o z -\- 



pb(s)ds. Next, we 



separate the spatial and temporal terms in velocity v(x,t) = 
v x (x)v t (t) and substitute the nano- and meso-stresses into the 
truncated constitutive equations (2) and (3). The nano-scale spatial 
perturbations v x (x) in the liquid and the solid are Gaussian 

and, then, become harmonic 



\x)=v z exp[-(x-z) 2 /2l 2 nan ^ 



Zij^CijklGkl 



(2) 



at meso-scale u™ eso {x) = 8v z Re[exp(i(x — z)/l meso )\, where So z is 
the magnitude of the perturbations. In the next step, we determine 
the temporal contribution v t (t). It turns out that the expression for 
v t (t) is scale independent (regardless whether it is the nano- or meso- 
scale), but v t (t) is different in the liquid and the solid 
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D^ md (t)=Re 
vf M {t)=Re 



' ' 2x 0 {Rf) 2 1 ' 



In here, the characteristic relaxation time is Tq = Vi c /v^, while 

v s — \/C/p is sound velocity and, as before, = rj/p is kinematic 
viscosity. The elastic constant is reduced to a single parameter C. A 
substance in the liquid and the solid state has very different 
properties and these properties become comparable only within 
the LST regime. The two expressions in (4) become identical when 
the material Reynolds number R^f is equal to one. We conclude that 
the liquid-solid transition emerges when 

° M = o,t/o lt = l. (5) 



R 1 



The magnitude of kinematic viscosity increases as the size of the 
confinement becomes smaller 7,837,38 . As indicated, the transitional 
length in water, hexadecane, cyclohexane and other substances is 
in the range of six to ten molecular distances. When knowing 
sound velocity and viscosity (both measured at small scales), the 
predicted characteristic length l c = i)k/v s is within the range 
observed in the experiments and predicted in molecular dynamics 
simulations. 

LLT- steady flow. We have shown that spatial nano-fluctuations in 
equilibrium are Gaussian and become harmonic at meso-scale. It is 
suggested that steady shear makes the fluctuations non-Gaussian 39 
and, then, the fluctuations (1) trigger the liquid-liquid transition. We 
begin by enforcing conservation of mass 8p/8t + 8(pv) /8x = 0 and 
momentum da/dx = pu + pudv/x, where the Cauchy stress is r/, 
velocity is u = du/dt and displacement is denoted as u. All the 
variables are expressed in terms of moving coordinate system 
g = x — Dt, where the steady velocity is D. Consequently, we have 
« = «(?)> a = (j(q) and strain is e = du/dq. Stress is derived directly 
from the conservation laws and is o(q) = — p 0 Do(cJ. In here, velocity 
is 0(c) = —Ddu/dg and p 0 is the initial mass density. 

As in the near equilibrium scenario, we predict the liquid-solid 
transition by comparing Eqns. 2 and 3; in here, both the relations 
include the contributions of the dynamic heterogeneities. The trans- 
ition occurs when 



#=£f- = l or - = 1. 



D 



D 



(6) 



Note that at D = D s the two material numbers are equal (R s e =Rf) 
and, consequently, the two criteria (5) and (6) become identical. The 
transitional liquid (liquid II) should exhibit the properties described 
by (2) and (3). From the solution presented in Methods A, flow 
patterns exhibited by the liquid II are limited to 



»(c) = 



TO 



"i«7M l)lti \ K - 1 



-u 2 e. 



RS-1 



where u x and u 2 are constants. From (7) we see that the liquid II 
cannot be formed at R s e = 1 . In all other situations (_R^ ^ 1 ) , the liquid 
II must follow the script defined in (7). There are four scenarios: 

(1) At R s e < 1, the substance is in its solid state (pastes, wet sands, 
and other dense suspensions) exhibiting high resistance to flow 
(high viscosity). As a result, the material Reynolds number is 
smaU (Rf < l) . The solid-liquid conversion is accomplished by 
applying a standing wave designed according to the protocol 



(3) 



(7). The most known phenomenon of this kind is soil 
liquefaction 40 . 

(2) One may envision another situation, where a viscoelastic liquid 
(R s e > l) is subjected to the standing wave. In this manner, the 
liquid is forced to act as if it were a liquid II substance. Such 
experiments have been conducted 14,41 and show the formation 
of microstructural patterns. We are interested in predicting the 
aggregation and disaggregation of red cells in blood not only 
under the steady state but also during transient flow 42 . 
Exponential flow of liquid II occurs at R^ > 1. Strongly driven 
liquids (under shock or impact) fit well the scenario. Often, it is 
assumed that shock wave (D > t) s ) has a sharp shock transition. 
In reality, the transition consists of several molecular layers of 
aligned molecules 15 which (we predict) are organized within a 
thin membrane. The membrane travels through the material 
with shock velocity. Impact loading is also known to trigger 
the liquid-liquid conversions 16 . 
(4) We predict that the liquid II behavior may emerge in liquids 
pushed from equilibrium and, then, allowed to relax according 
to the rules in (7). Supercooled liquids may fit the scenario, 
where a controllable decrease of temperature leads to an 
increase of viscosity, thus, affecting the substance's relaxation 
time 21 . We are not aware of any other experimental work done 
in this area. 

LLST- short-lived transitions in blood during cardiac cycle. In 

simple terms, blood is a liquid tissue consisting of plasma and 
blood cells. On average 1 microliter of blood contains about 5T0 6 
red cells. Thus, in larger vessels (diameter 2 mm or larger), the 
number of cells is large enough for blood to become a homo- 
geneous viscoelastic liquid. Blood viscosity and elasticity strongly 
depend on the actual blood composition, flow rate, shape and size 
of the blood vessels. A steady flow is fastest at the center of the vessel 
and slowest near the wall. This non- uniformity is linked to the 
buildup of wall shear stress. A pulsatile cycle produces conside- 
rably more complex flow patterns 43,44 . It is observed that vessel 
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Figure 2 | A rigid tube filled with blood is rapidly moved in axial direction 
(x) at constant velocity. Initially, the velocity gradient triggers liquid- 
liquid-solid transitions but, later, blood is converted back to its liquid 
form. The liquid (L), transitional liquid (T) and solid (S) are marked in the 
tube and are plotted as a function of time. Velocity distributions are shown 
in 2-D and 3-D representations and are plotted for t = 0, 50, 100 and 
150 ms. 
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Figure 3 | In the upper right corner of the figure, we show the relaxation of the solid (S) and transitional (T) phases as a function of time for three vessels 
(D = 3.78, 4.54 and 6.04 nm). The liquid II (T) and the solid (S) phases vanish faster as the size of the vessel increases. Also, the thickness of the T-S layer is 
plotted as a function of the vessel diameter. At the initial time of the flow reversal ( t = 0) , the T-S layer has a measurable thickness even in a very large vessel 
(D > 20 mm). However, the size of the layer shrinks with the progression of time, as shown in instances (t = 10, 20, 30, 40 and 50 ms). 



segments with low wall shear stress and oscillatory changes in flow 
direction appear to be at high risk for the development of various 
diseases and among them atherosclerosis and thrombosis. 
Atherosclerosis affects the inner lining of an artery and is 
characterized by plaque deposits that block the flow of blood. 
Thrombosis is the formation of clots capable of obstructing the 
flow of blood. Factors that affect blood viscosity 45 " 48 are: 
hematocrit, red cell aggregation and deformation, plasma viscosity, 
concentration and size of low-density lipoproteins (LDL) and age. 

Our objective is to forecast liquid-liquid-solid transitions in blood 
during cardiac cycle 49 . In our idealized scenario, the vessel is a rigid 
tube with inner radius l e . The tube is filled with blood and the system 
stays at rest. The derivations are based on Lagrange description, 
where strains are small, times are relatively short t < t 0 . The tube is 
large enough for the blood to remain in its liquid state (Rf>2). 
Next, the tube is rapidly moved from rest along its axis and kept in 
motion at constant velocity Vo, Fig. 2. Blood slippages along the walls 
are not allowed. Whole blood in the state (L) is the viscoelastic liquid 
(2), where the temporal and spatial fluctuations are considered 
important. The blood in the transitional state (T) exhibits the prop- 
erties of both the liquid (L) and the solid (S). The solid-like behavior 
is described in terms of the "local" Kelvin- Voigt model, where the 
stress fluctuations are omitted. The liquid (L), liquid II (T) and solid 
(S) are glued together, Fig. 2. The boundary of the liquid r L with 
respect to the boundary of the solid r s in the presence of the liquid 
II ( rs — Tl ) must be optimal in terms of the rate of work performed by 
the system (owv 0 ), where the wall shear stress is aw = a(l e ,t) and 
a(r,t ) is shear stress. Expressions for the particle velocity and stress in 
each state are presented in Methods B. Blood viscosity and elasticity 
are determined for whole-blood (hematocrit 38%), where 
v s — A. 3-10 _ 3 mmjms and Vk = 3. 25-10 ~ 3 mm 2 /ms. In each solution, 
the radius l e of the vessel is equal to the significant stretch l c of the 
stress gradient. As stated earlier, the vessel is rapidly moved from rest 
and kept in motion at constant velocity t> 0 = 0.12mm/ms. In terms of 
cardiac cycle this is the worst case scenario. The analysis is con- 
structed for relatively large vessels, where the diameter is varying 
between three to six millimeters. These diameters correspond to 
the material Reynolds number R 1 ^ in the range of two to four. 
Rapid departure from rest converts the liquid (L) into the liquid II 
(T) and the solid-like material (S), Fig. 2. In all the studied cases, clock 
is set to zero when all particles across the tube start sensing the 



motion, while the particle velocity at the center of the tube is still 
equal to zero. This setup properly replicates the velocity distribution 
in the tube at the moment of the blood flow reversal 49 . In the first 
example (i?^ = 2) , at r = 0 the entire vessel contains blood in the 
transitional (T) and solid (S) states, Fig. 2. Gradually, blood is con- 
verted back to its liquid form (L) and at t ~ z 0 the conversion is 
complete, where t 0 = 176ms. There is a moment when the trans- 
itional liquid (T) disappears and a sharp liquid-solid interface 
emerges (rs = ri). At this point, the L-S interface migrates toward 
the tube walls. Stress tractions along the L-S interface are satisfied but 
velocity becomes discontinuous causing slippages. Such slippages 
have been observed in vessels near the vessel walls in a plasma layer 
of the thickness about 45.8 p:m 15 ' 49 . In larger vessels(^ = 2.5,3,4), 
the layers T and S are smaller and the conversion process is faster 
(Fig. 3). An indirect support for the LLSTs in blood at transient 
conditions is offered in ref. 15. The presence of the transitional 
and solid-like blood near the vessel walls is a concerning factor. 
We should note that an increase of blood viscosity and/or measurable 
decrease of elasticity may further aggravate the problem. Often, 
blood viscosity is considered the unifying indicator of cardiovascular 
diseases. It seems that the material Reynolds number would be a 
better predictor of cardiovascular disease risk. 

Discussion 

There are two aspects of the work worth noticing. First, our con- 
tinuum level approach is adapted for a small scale analysis. We 
accomplish this by incorporating the spatial and temporal stress 
fluctuations into the material's description. Second, we view the 
substance either as a liquid-like or solid-like material. Thus, we 
diverge from the approaches where the liquid-like, solid-like and 
the transitional behaviors are constructed within a single mathemat- 
ical framework. In this manner, we gain additional degree of freedom 
which we find necessary when describing the transitional processes. 
Consequently, we predict LLSTs during active flow processes (steady 
or transient); where in some circumstances the substance is hard 
driven. There are several mechanisms by which the transitions occur. 
Our criterion for the liquid-solid transformation works for liquids 
in small confinements. We predict microstructural reorganiza- 
tions in liquids stimulated by standing waves, as we show that the 
standing waves are responsible for triggering liquefaction in dense 
suspensions. A synchronized motion of molecules (or particles in 
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suspension) is forecasted in liquids subjected to an impact and/or 
shock loading. Lastly, we indicate that the liquid-liquid transitions 
should occur in liquids pushed away from equilibrium and, then, are 
allowed to relax according to a controllable scenario. 

Methods 

A. Transitional liquid (liquid II): 

Suppose that the liquid-like and solid-like behaviors are indistinguishable. If 
so, we can substitute stress from (3) into the flow equation (2) and construct a general 
solution. First, we introduce an effective velocity 

u(x,t) du(x,t) 



~v(x,t) -- 



To 



dt 



Next, we define an effective strain rate 

du(x,t) t 0 db{x,t) 



s(x,t) = - 



i /.V 



dl 



(A.l) 



(A.2) 



The liquid II obeys the relations (2) and (3). We find that with the use of the new 
variables (A.l) and (A.2) the transitional liquid must follow the rule defined below 



e(*,f) + h — - y [o(x,t 0 ) + l c £(x,to)]. 
ax l c 



(A3) 



In here, o(x s to) is the initial velocity and e(x4o) is the initial strain rate. In addition, 
the material in the transitional state must adhere to the conservation laws. Under the 
condition of steady flow (q — x — Dt), the rules (A3) are 



d 

3c 



I"0 TO 



To \ du(q) 
k J dc 



ID 



d 2 u(g) 
d? 



-D 



t(^ +tD d Q - 



d 2 u(g) 



(A.4) 



where stress is cr(cj — p Q D 2 -du(g) / dq. Solution of the problem is presented in terms of 
displacements and stresses, namely 



u(q) = u a + u l exp\^-J \] RS-lj +u 2 exp 
a M = Dp ^_ L exp (-5- 1 Jr?-J 



_ C 

U2CXp 



(A.5) 



B. LLST in rigid vessel: 

In a rigid vessel filled with blood, velocity and shear stress in the liquid (L) 
satisfy the equation of motion and are 

2Rf(Rfr + 2o s t) 



o L ( r t) = o a , +«) I - H — 

\2 (R») 



fl ex P 



1 1 P D * I , 



[St 



■j\exp[ 



2Kf (fif r+2v s t) 



(B.l) 



(4+(«f) 2 



There are two constants, namely o° L and o l L . Moreover, the velocity gradient in the 
center of the tube is always equal to zero Sl(O^) —0, where s L — du L /dr. 

The responses of the liquid II are constructed by solving the equation (A.3). 
Consequently, velocity and shear stress are 

VT (r,t) = v 0 T + u l T fexp[-(r + D s t)/l e } + i< 2 T ^-^exp[-R™ (R?r+o s t) /l e ] 

(B.2) 

o T = - { u\exp[- (r + v s t) /Q + ii^exp [ - Rf (i?f r + a,t) //J } ^ . 

The transitional flow is determined in terms of three constants u° , u l T and u\. 
Lastly, the solution for the solid is 

D S (r,t) = u° + 4r 

I \ .„ (B-3) 



where u° s and iP s are constants. Boundary conditions for this problem are defined as 
follows: 



o s (U) = o„; o L ft,0) = o<,; »l(0,0) =0; 
Os(rs,t)=i>r(rs,ty,<rs(rs,t)=<TT(rs,ty, 
o L (r L ,t) = o T (r L ,ty,a L (r L ,t) = a T (r L ,t). 



(B.4) 



We have seven constants {(;^,e°,(;°,tj},(j^,w^,u^ j and two time-dependent variables 
{ r L( f )' r s(0}- The LLT and LLST boundaries are determined from the criterion of 
least action 



d[<Tw{rL-,r$,tyvo]/dr s = 0. 



(B.5) 
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